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Abstract 

We describe a major update of our Matlab freeware GloptiPoly for parsing gen- 
eralized problems of moments and solving them numerically with semidefinite pro- 
gramming. 



1 What is GloptiPoly ? 



> 

On , 

, Gloptipoly 3 is intended to solve, or at least approximate, the Generalized Problem of 

Moments (GPM), an infinite-dimensional optimization problem which can be viewed as 

ON • 

an extension of the classical problem of moments [8]. From a theoretical viewpoint, the 

: 

GPM has developments and impact in various areas of mathematics such as algebra, 
Fourier analysis, functional analysis, operator theory, probability and statistics, to cite 
a few. In addition, and despite a rather simple and short formulation, the GPM has a 
large number of important applications in various fields such as optimization, probability, 
finance, control, signal processing, chemistry, cristallography, tomography, etc. For an 
account of various methodologies as well as some of potential applications, the interested 
reader is referred to [lj [2] and the nice collection of papers [5]. 

The present version of GloptiPoly 3 can handle moment problems with polynomial data. 
Many important applications in e.g. optimization, probability, financial economics and 
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optimal control, can be viewed as particular instances of the GPM, and (possibly after 
some transformation) of the GPM with polynomial data. 

The approach is similar to that used in the former version 2 of GloptiPoly Sj. The 
software allows to build up a hierarchy of semidefmite programming (SDP), or linear 
matrix inequality (LMI) relaxations of the GPM, whose associated monotone sequence of 
optimal values converges to the global optimum. For more details on the approach, the 
interested reader is referred to [8]. 

2 Installation 

GloptiPoly 3 is a freeware subject to the General Public Licence (GPL) policy. It can be 
downloaded at 

www. laas . f r/~henrion/ sof tware/gloptipoly3 

The package, available as a compressed archive, consists of several m-files and subdirec- 
tories, and it contains no binaries. Extracted files are placed in a gloptipoly3 directory 
that should be declared in the Matlab working path, using e.g. Matlab's command 

» addpath gloptipoly3 

GloptiPoly 3 uses by default the semidefmite programming solver SeDuMi [TTj . so this 
package should be properly installed. Other semidefinite solvers can also be used provided 
they are installed and interfaced through YALMIP |10j . 

3 Getting started 

Please type the command 
» gloptipolydemo 
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to run interactively the basic example that follows. 

Consider the classical problem of minimizing globally the two-dimensional six-hump camel 
back function [3] 

min go(x) = Ax\ + x\Xi — Ax\ — 2.1x\ + Ax\ + -x\. 

xeu 2 3 

The function has six local minima, two of them being global minima. 

Using GloptiPoly 3, this optimization problem can be modeled as a moment problem as 
follows: 

» mpol xl x2 

» g = 4*xl~2+xl*x2-4*x2~2-2. l*xl~4+4*x2~4+xl~6/3 
Scalar polynomial 

4xl~2+xlx2-4x2~2-2 . lxl"4+4x2"4+0 . 33333xl~6 

» P = msdp(min(gO)) ; 

GloptiPoly 3.0 

Define moment SDP problem 

(GloptiPoly output suppressed) 

Generate moment SDP problem 
» P 

Moment SDP problem 

Measure label = 1 

Relaxation order = 3 

Decision variables = 27 

Semidefinite inequalities = 10x10 

Once the moment problem is modeled, a semidefinite solver can be used to solve it nu- 
merically. Here we use SeDuMi [11] which is assumed to be installed and accessible from 
the Matlab working path: 
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» [status, obj] = msol(P) 

GloptiPoly 3.0 

Solve moment SDP problem 

Calling SeDuMi 

SeDuMi 1.1R3 by AdvOL , 2006 and Jos F. Sturm, 1998-2003. 

(SeDuMi output suppressed) 

2 globally optimal solutions extracted 
» status 
status = 
1 

» obj 
obj = 

-1.0316 
» x = double ( [xl x2] ) ; 
x(:,:,l) = 

0.0898 -0.7127 
x(:,:,2) = 

-0.0898 0.7127 

The flag status = 1 means that the moment problem is solved successfully and that 
GloptiPoly can extract two globally optimal solutions reaching the objective function obj 
= -1.0316. 

4 From version 2 to version 3 

The major changes incorporated into GloptiPoly when passing from version 2 to 3 can be 
summarized as follows: 
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• Use of native polynomial objects and object-oriented programming with specific 
classes for multivariate polynomials, measures, moments, and corresponding over- 
loaded operators. In contrast with version 2, the Symbolic Toolbox for Matlab 
(gateway to the Maple kernel) is not required anymore to process polynomial data. 

• Generalized problems of moments featuring several measures with semialgebraic 
support constraints and linear moment constraints can be processed and solved. 
Version 2 was limited to moment problems on a unique measure without moment 
constraints. 

• Explicit moment substitutions are carried out to reduce the number of variables and 
constraints. 

• The moment problems can be solved numerically with any semidefinite solver, pro- 
vided it is interfaced through YALMIP. In contrast, version 2 used only the solver 
SeDuMi. 

5 Solving generalized problems of moments 

GloptiPoly 3 uses advanced Matlab features for object-oriented programming and over- 
loaded operators. The user should be familiar with the following basic objects. 

5.1 Multivariate polynomials (mpol) 

A multivariate polynomial is an affine combination of monomials, each monomial de- 
pending on a set of variables. Variables can be declared in the Matlab working space as 
follows: 

» clear 
» mpol x 
» x 

Scalar polynomial 
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X 

» mpol y 2 

» y 

2- by-l polynomial vector 

(i,D:y(D 

(2,1) :y(2) 
» mpol z 3 2 
» z 

3- by-2 polynomial matrix 
(1,1) :z(l,l) 

(2,1) :z(2,l) 

(3.1) :z(3,l) 

(1.2) :z(l,2) 
(2,2) :z(2,2) 
(3,2) :z(3,2) 

Variables, monomials and polynomials are defined as objects of class mpol. 
All standard Matlab operators have been overloaded for mpol objects: 

» y*y ' -z ' *Z+X~3 

2-by-2 polynomial matrix 

(1,1) :y(l)~2-z(l,l)-2-z(2,l)~2-z(3,l)-2+x~3 

(2.1) :y(l)y(2)-z(l,l)z(l,2)-z(2,l)z(2,2)-z(3,l)z(3,2)+x~3 

(1.2) :y(l)y(2)-z(l,l)z(l,2)-z(2,l)z(2,2)-z(3,l)z(3,2)+x~3 
(2,2) :y(2)"2-z(l,2)"2-z(2,2)"2-z(3,2)"2+x~3 

Use the instruction 

» mset clear 

to delete all existing GloptiPoly variables from the Matlab working space. 
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5.2 Measures (meas) 



Variables can be associated with real-valued measures, and one variable is associated with 
only one measure. For GloptiPoly, measures are identified with a label, a positive integer. 
When starting a GloptiPoly session, the default measure has label 1. By default, all 
created variables are associated with the current measure. Measures can be handled with 
the class meas as follows: 

» mset clear 
» mpol x 
» mpol y 2 
» meas 

Measure 1 on 3 variables: x,y(l),y(2) 
» meas(y) °/ create new measure 
Measure 2 on 2 variables: y(l),y(2) 
» m = meas 

l-by-2 vector of measures 
1: Measure 1 on 1 variable: x 
2:Measure 2 on 2 variables: y(l),y(2) 
» m(l) 

Measure number 1 on 1 variable: x 

The above script creates a measure dfj,i(x) on M and a measure d^iy) on IR 2 . 
Use the instruction 

» mset clearmeas 

to delete all existing GloptiPoly measures from the working space. Note that this does 
not delete existing GloptiPoly variables. 
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5.3 Moments (mom) 



Linear combinations of moments of a given measure can be manipulated with the mom 
class as follows: 

» mom(l+2*x+3*x~2) 
Scalar moment 
I[l+2x+3x~2]d[l] 
» mom(y*y ) ) 
2-by-2 moment matrix 
(1,1) :I[y(l)~2]d[2] 

(2.1) :I[y(l)y(2)]d[2] 

(1.2) :I[y(l)y(2)]d[2] 
(2,2) :I[y(2)~2]d[2] 

The notation I [p] d [k] stands for f p where p is a polynomial of the variables asso- 
ciated with measure d/ik, and k is the measure label. 

Note that it makes no sense to define moments over several measures, or nonlinear moment 
expressions: 

» mom(x*y(l)) 

??? Error using ==> mom. mom 

Invalid partitioning of measures in moments 

>> mom(x) *mom(y(l) ) 

??? Error using ==> mom. times 

Invalid moment product 

Note also the distinction between a constant term and the mass of a measure: 



» l+mom(x) 
Scalar moment 
1+I[x]d[l] 
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» mom(l+x) 
Scalar moment 
I[l+x]d[l] 
» mass(x) 
Scalar moment 
I[l]d[l] 



Finally, let us mention three equivalent notations to refer to the mass of a measure: 

» mass(meas(y)) 
Scalar moment 
I[l]d[2] 
» mass(y) 
Scalar moment 
I[l]d[2] 
» mass (2) 
Scalar moment 
I[l]d[2] 



The first command refers explicitly to the measure, the second command is a handy short- 
cut to refer to a measure via its variables, and the third command refers to GloptiPoly's 
labeling of measures. 



5.4 Support constraints (supcon) 

By default, a measure on n variables is defined on the whole M™. We can restrict the 
support of a mesure to a given semialgebraic set as follows: 

» 2*x~2+x~3 == 2+x 

Scalar measure support equality 

2x~2+x~3 == 2+x 

» disk = (y'*y <= 1) 
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Scalar measure support inequality 
y(l)~2+y(2)~2 <= 1 



Support constraints are modeled by objects of class supcon. The first command means 
that variable x must satisfy x 3 + 2x 2 — x — 2 = (x — l)(x + l)(x + 2) =0, i.e. measure 
dfj,i(x) must be discrete, a linear combination of three Dirac at 1, —1 and —2. The second 
command restricts measure dfj, 2 (y) within the unit disk. 

Note that it makes no sense to define a support constraint on several measures: 
» x+y(l) <= 1 

??? Error using ==> supcon. supcon 
Invalid reference to several measures 

5.5 Moment constraints (momcon) 

We can constrain linearly the moments of several measures: 

» mom(x~2+2) == l+mom(y (1) "3*y (2) ) 
Scalar moment equality constraint 
I[2+x~2]d[l] == 1+I[y(l)~3y(2)]d[2] 
» mass(x)+mass(y) <= 2 
Scalar moment inequality constraint 
I[l]d[l]+I[l]d[2] <= 2 

Moment constraints are modeled by objects of class momcon. 

For GloptiPoly an objective function to be minimized or maximized is considered as a 
particular moment constraint: 

» min(mom(x~2+2)) 

Scalar moment objective function 

min I[2+x~2]d[l] 

10 



» max(x"2+2) 

Scalar moment objective function 
max I[2+x~2]d[l] 



The latter syntax is a handy short-cut which directly converts an mpol object into an 
momcon object. 

5.6 Floating point numbers (double) 

Variables in a measure can be assigned numerical values: 

» ml = assign(x,2) 
Measure 1 on 1 variable: x 
supported on 1 point 

which is equivalent to enforcing a discrete support for the measure. Here d\x\ is set to the 
Dirac at the point 2. 

The double operator converts a measure or its variables into a floating point number: 

» double (x) 
ans = 

2 

» double (ml) 
ans = 

2 

Polynomials can be evaluated similarly: 

»double ( l-2*x+3*x~2) 
ans = 

9 
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Discrete measure supports consisting of several points can be specified in an array: 

» m2 = assign(y, [-1 2 0;l/3 1/4 -2]) 
Measure 2 on 2 variables: y(l) ,y(2) 
supported on 3 points 
» double (m2) 
ans( : , : , 1) = 
-1.0000 

0.3333 
ans( : , : ,2) = 

2.0000 

0.2500 
ans( : , : ,3) = 


-2 

5.7 Moment SDP problems (msdp) 

GloptiPoly 3 can manipulate and solve Generalized Problems of Moments (GPM) as 
defined in [8]: 



In the above notations, gik{x), hjk(x) are given real polynomials and bj are given real 
constants. The decision variables in the GPM are measures d/J,k(x), and GloptiPoly 3 
allows to optimize over them through their moments 




where measures d\x^ are supported on basic semialgebraic sets 



K fc = {^R"' = 9ik(x) >0, i = l,2...}. 




where the are multi-indices. 
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5.8 Solving moment problems msol 



Once a moment problem is defined, it can be solved numerically with the instruction 
msol. In the sequel we give several examples of GPMs handled with GloptiPoly 3. 

5.8.1 Unconstrained minimization 

Following [6], given a multivariate polynomial go(x), the unconstrained optimization prob- 
lem 

min qo(x) 

can be formulated as a linear moment optimization problem 

min^ / g (x)d/i(x) 
s.t. / dfi(x) = 1 

where measure d\i lives in the space B n of finite Borel signed measures on M n . The 
equality constraint indicates that the mass of d\x is equal to one, or equivalently, that d/i 
is a probability measure. 

In general, this linear (hence convex) reformulation of a (typically nonconvex) polynomial 
problem is not helpful because there is no computationally efficient way to represent 
measures and their underlying Borel spaces. The approach proposed in [6] consists in using 
convex semidefinite representations of the space B n truncated to finite degree moments. 
GloptiPoly 3 allows to input such moment optimization problems in an user-friendly way, 
and to solve them using existing software for semidefinite programming (SDP). 

In Section [3] we already encountered an example of an unconstrained polynomial opti- 
mization solved with GloptiPoly 3. Let us revisit this example: 

>> mset clear 
» mpol xl x2 

» g = 4*xl~2+xl*x2-4*x2~2-2. l*xl~4+4*x2~4+xl~6/3 
Scalar polynomial 

4xl~2+xlx2-4x2~2-2 . lxl"4+4x2"4+0 . 33333xl~6 
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» P = msdp(min(gO)) ; 
» msol(P) 

2 globally optimal solutions extracted 
Global optimality certified numerically 

This indicates that the global minimum is attained with a discrete measure supported on 
two points. The measure can be constructed from the knowledge of its first moments of 
degree up to 6: 

» meas 

Measure 1 on 2 variables: xl,x2 

with moments of degree up to 6, supported on 2 points 
» double (meas) 
ans( : , : , 1) = 
0.0898 
-0.7127 
ans( : , : ,2) = 
-0.0898 
0.7127 
» double (gO) 
ans( : , : , 1) = 

-1.0316 
ans( : , : ,2) = 
-1.0316 

When converting to floating point numbers with the operator double, it is essential to 
make the distinction between mpol and mom objects: 

» v = mmon( [xl x2] ,2) ' 
l-by-6 polynomial vector 
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(1.1) :1 

(1.2) :xl 

(1.3) :x2 

(1.4) :xl~2 

(1.5) :xlx2 

(1.6) :x2"2 
» double (v) 
ans( : , : , 1) = 

1.0000 0.0898 -0.7127 0.0081 -0.0640 0.5079 
ans( : , : ,2) = 

1.0000 -0.0898 0.7127 0.0081 -0.0640 0.5079 
» double (mom (v)) 
ans = 

1.0000 0.0000 -0.0000 0.0081 -0.0640 0.5079 

The first instruction mmon generates a vector of monomials v of class mpol, so the command 
double (v) calls the convertor @mpol/double which evaluates a polynomial expression on 
the discrete support of a measure (here two points). The last command double (mom (v)) 
calls the convertor @mom/double which returns the value of the moments obtained after 
solving the moment problem. 

Note that when inputing moment problems on a unique measure whose mass is not con- 
strained, GloptiPoly assumes by default that the measure has mass one, i.e. that we are 
seeking a probability measure. Therefore, if gO is the polynomial defined previously, the 
two instructions 

» P = msdp(min(g0)) ; 

and 

» P = msdp(min(g0) , mass(meas(g0))==l) ; 

are equivalent. See also Section l5\3l for handling masses of measures and Section [5. 8. 21 for 
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more information on mass constraints. 



5.8.2 Constrained minimization 



Consider now the constrained polynomial optimization problem 



min g (x) 



where 



K = {x G R n : gi(x) > 0, i = 1, 2, . . . 



} 



is a basic semialgebraic set described by given polynomials g%{x). Following [6], this (non- 
convex polynomial) problem can be formulated as the (convex linear) moment problem 



where the indeterminate is a probability measure d\i of B n which is now supported on set 
K. In other words 



JR"/K 

As an example, consider the non-convex quadratic problem of Section 4.4 in [3]: 
min — 2x\ + x 2 — x 3 

s.t. 24 - 20xi + 9x 2 - 13x 3 + Ax{ - Ax x x 2 + 4xix 3 + 2x\ - 2x 2 x 3 + 2x\ > 
x 1 + x 2 + x 3 < 4, 3x 2 + x 3 < 6 
< xi < 2, < x 2 , < x 3 < 3 

Each constraint in this problem is interpreted by GloptiPoly 3 as a support constraint on 
the measure associated with variable x, see Section 15.41 

» mpol x 3 

» x(l)+x(2)+x(3) <= 4 

Scalar measure support inequality 

x(l)+x(2)+x(3) <= 4 

The whole problem can be entered as follows: 



min dM J K g (x)d^(x) 
s.t. j K dn(x) = l 
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» mpol x 3 

» g o = -2*x(l)+x(2)-x(3); 

»K= [24-20*x(l)+9*x(2)-13*x(3)+4*x(l)"2-4*x(l)*x(2) ... 

+4*x(l)*x(3)+2*x(2)~2-2*x(2)* x (3)+2*x(3)~2 >= 0, ... 
x(l)+x(2)+x(3) <= 4, 3*x(2)+x(3) <= 6, ... 
<= x(l), x(l) <= 2, <= x(2), <= x(3), x(3) <= 3]; 
» P = msdp(min(gO) , K) 

Moment SDP problem 

Measure label = 1 

Relaxation order = 1 

Decision variables = 9 

Linear inequalities = 8 

Semidefinite inequalities = 4x4 

The moment problem can then be solved: 

» [status, obj] = msol(P) 

GloptiPoly 3.0 

Solve moment SDP problem 

Global optimality cannot be ensured 
status = 


obj = 

-6.0000 

Since status=0 the moment SDP problem can be solved but it is impossible to detect 
global optimality. The value obj =-6. 0000 is then a lower bound on the global minimum 
of the quadratic problem. 

The measure associated with the problem variables can be retrieved as follows: 
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» mu = meas 

Measure 1 on 3 variables: x(l) ,x(2) ,x(3) 
with moments of degree up to 2 

Its vector of moments can be built as follows: 

» mv = mvec(mu) 
10-by-l moment vector 
(l,l):I[l]d[l] 
(2,l):I[x(l)]d[l] 
(3,l):I[x(2)]d[l] 
(4,l):I[x(3)]d[l] 
(5,1) :I[x(l)~2]d[l] 
(6,1) :I[x(l)x(2)]d[l] 
(7,l):I[x(l)x(3)]d[l] 
(8,1) :I[x(2)~2]d[l] 
(9,1) :I[x(2)x(3)]d[l] 
(10,1) :I[x(3)"2]d[l] 

These moments are the decision variables of the SDP problem solved with the above msol 
command. Their numerical values can be retrieved as follows: 

» double (mv) 
ans = 

1.0000 

2.0000 
-0.0000 

2.0000 

7.6106 

1.4671 

2.3363 

4.8335 
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0.5008 
8 . 7247 



The numerical moment matrix can be obtained using the following commands: 

» double (mmat(mu)) 
ans = 



1. 


.0000 


2. 


.0000 


-0. 


.0000 


2 


.0000 


2. 


.0000 


7. 


.6106 


1. 


.4671 


2 


.3363 


-0. 


,0000 


1. 


.4671 


4. 


.8335 





.5008 


2. 


.0000 


2. 


.3363 


0. 


.5008 


8 


.7247 



As explained in |6j, we can build a hierarchy of nested moment SDP problems, or relax- 
ations, whose solutions converge monotically and asymptotically to the global optimum, 
under mild technical assumptions. By default the command msdp builds the relaxation of 
lowest order, equal to half the degree of the highest degree monomial in the polynomial 
data. An additional input argument can be specified to build higher order relaxations: 

» P = msdp(min(g0) , K, 2) 

Moment SDP problem 

Measure label = 1 

Relaxation order = 2 

Decision variables = 34 

Semidefinite inequalities = 10xl0+8x(4x4) 
» [status, obj] = msol(P) 

Global optimality cannot be ensured 
status = 







obj = 



-5.6922 
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» P = msdp(min(gO) , K, 3) 

Moment SDP problem 

Measure label = 1 

Relaxation order = 3 

Decision variables = 83 

Semidefinite inequalities = 20x20+8x(10xl0) 
» [status, obj] = msol(P) 

Global optimality cannot be ensured 
status = 


obj = 

-4.0684 

We observe that the moment SDP problems feature an increasing number of variables 
and constraints. They generate a mononotically increasing sequence of lower bounds on 
the global optimum, which is eventually reached numerically at the fourth relaxation: 

» P = msdp(min(gO) , K, 4) 

Moment SDP problem 
Measure label 
Relaxation order 
Decision variables 
Semidefinite inequalities 
» [status, obj] = msol(P) 

2 globally optimal solutions extracted 
Global optimality certified numerically 
status = 
1 



= 1 
= 4 
= 164 

= 35x35+8x (20x20) 
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obj = 

-4.0000 

» double (x) 

ans( : , : , 1) = 
2.0000 
0.0000 
0.0000 

ans( : , : ,2) = 
0.5000 
0.0000 
3.0000 

» double (gO) 

ans( : , : , 1) = 
-4.0000 

ans( : , : ,2) = 
-4.0000 



5.8.3 Rational minimization 

Minimization of a rational function can also be formulated as a linear moment problem. 

Given two polynomials go{x) and ho(x), consider the rational optimization problem 

. 9o(x) 
mm 

where 

X = {x G M n : gi(x) > 0, i = 1, 2, . . .} 

is a basic semialgebraic set described by given polynomials gi(x). Following [4J, the 
corresponding moment problem is given by 

min dMe B™ f K g {x)dfJ,(x) 

s.t. L ho(x)dfj,(x) = 1. 

In contrast with the polynomial optimization problem of Section 15.8.21 the optimal mea- 
sure dfi supported on X is not necessarily a probability measure. Denoting h (x) = 
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Yla hoaX a , the moments y a of d\i must satisfy a linear constraint 

/ ho(z)d(i(x) = y2h 0a y a = 1. 

As an example, consider the one- variable rational minimization problem [H Ex. 2]: 



mm 



+ 2x + 1 

We can solve this problem with GloptiPoly 3 as follows: 
» mpol x 

» gO = x~2-2*x; hO = x~2+2*x+l; 
» P = msdp(min(gO) , mom(hO) == 1); 
» [status, obj] = msol(P) 

Global optimality certified numerically 
status = 
1 

obj = 

-0.3333 
» double (x) 
ans = 

0.4999 
5.8.4 Several measures 

GloptiPoly 3 can handle several measures whose moments are linearly related. 

For example, consider the GPM arising when solving polynomial optimal control problems 
as detailed in [7]. We are seeking two occupation measures dfj,i(x, u) and d^ix) of a state 
vector x(t) and input vector u(t) whose time variation are governed by the differential 
equation 

^JT = f(x,u), x(0) = x , u(0) = m 
dt 
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with f(x, u) a given polynomial mapping and x , u given initial conditions. Measure dfj,% 
is supported on a given semialgebraic set K x corresponding to constraints on x and u. 
Measure d^2 is supported on a given semialgebraic set K2 corresponding to performance 
requirements. For example K2 = indicates that state x must reach the origin. 

Given a polynomial test function g(x) we can relax the dynamics constraint with the 
moment constraint 

I g{x)dii 2 {x) - g(x ) = f ^ 1 X f(x,u)dni(x,u) 
Jk 2 Jki ''■>' 

linking linearly moments of dfii and d^2- As explained in |7|, a lower bound on the 
minimum time achievable by any feedback control law u(x) is then obtained by minimizing 
the mass of dfi\ over all possible measures d/j,\, d/j.2 satisfying the support and moment 
constraints. The gap between the lower bound and the exact minimum time is narrowed 
by enlarging the class of test functions g. 

In the following script we solve this moment problem in the case of a double integrator 
with state and input constraints: 

°/„ bounds on minimal achievable time for optimal control of 
% double integrator with state and input constraints 



xO = [1; 1]; uO = 0; % initial conditions 
d = 6; % maximum degree of test function 

% analytic minimum time 
if x0(l) >= -(x0(2)~2-2)/2 

tmin = 1+x0(1)+x0(2)+x0(2)-2/2; 
elseif x0(l) >= -x0(2)~2/2*sign(x0(2)) 

tmin = 2*sqrt(x0(l)+x0(2)~2/2)+x0(2) ; 
else 

tmin = 2*sqrt(-x0(l)+x0(2)~2/2)-x0(2) ; 
end 
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% occupation measure for constraints 
mpol xl 2 
mpol ul 

ml = meas( [xl ;ul] ) ; 

% occupation measure for performance 

mpol x2 2 

m2 = meas (x2) ; 

% dynamics 

scaling = tmin; % time scaling 
f = scaling* [xl (2) ;ul] ; 

°/„ test function 
gl = mmon(xl,d) ; 
g2 = mmon(x2,d) ; 

°/ initial condition 
assign ( [xl ;ul] , [xO;uO] ) ; 
gO = double (gl) ; 

% moment problem 

P = msdp(min(mass (ml) ) , . . . 
ul~2 <= 1,... % input constraint 
xl(2) >= -1,... °/„ state constraint 
x2'*x2 <= 0, . . . % performance = reach the origin 
mom(g2) - gO == mom(dif f (gl ,xl) *f ) ) ; °/ linear moment constraints 

% solve 

[status, obj] = msol(P); 
obj = scaling*obj ; 
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disp( ['Minimum time = ' num2str (tmin)] ) ; 

disp(['LMI ' int2str(d) ' lower bound = ' num2str (obj )] ) 

For the initial condition xq = [1 1] the exact minimum time is equal to 3.5. In Tabled] we 
report the monotically increasing sequence of lower bounds obtained by solving moment 
problems with test functions of increasing degrees. We used the above script and the 
semidefmite solver SeDuMi 1.1R3. 



degree 


2 


4 


6 8 


10 


12 


14 


16 


bound 


1.0019 


2.3700 


2.5640 2.9941 


3.3635 


3.4813 


3.4964 


3.4991 



Table 1: Minimum time optimal control for double integrator with state and input con- 
straints: lower bounds on exact minimal time 3.5 achieved by solving moment problems 
with test functions of increasing degrees. 

5.9 Using YALMIP 

By default GloptiPoly 3 uses the semidefmite solver SeDuMi [UJ for solving numeri- 
cally SDP moment problems. It is however possible to use any solver interfaced through 
YALMIP [10] by setting a configuration flag with the mset command: 

» mset (' yalmip true) 

Parameters for YALMIP, handled with the YALMIP command sdpsettings, can be 
forwarded to GloptiPoly 3 with the mset command. For example, the following com- 
mand tells YALMIP to use the SDPT3 solver (instead of SeDuMi) when solving moment 
problems with GloptiPoly: 

» mset (sdpsettings ( ' solver ' , ' sdpt3 ' ) ) ; 
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5.10 SeDuMi parameters settings 

The default parameters settings of SeDuMi [11] can be altered as follows: 

» pars.eps = le-10; 
» mset(pars) 

where pars is a structure of parameters consistent with SeDuMi's format. 

5.11 Exporting moment SDP problems 

A moment problem P of class msdp can be converted into SeDuMi's input format: 
» [A,b,c,K] = msedumi(P); 

The SDP problem can then be solved with SeDuMi as follows: 
» [x,y,info] = sedumi(A,b,c,K) ; 

See [11] for more information on SeDuMi's input data format. 

Similarly, a moment SDP problem can be converted into YALMIP's input format: 

» [F,h,y] = myalmip(P) ; 

where variable F contains the LMI constraints (YALMIP class lmi), h is the objective 
function (YALMIP class sdpvar) and y is the vector of moments (YALMIP class sdpvar). 
The SDP problem can then be solved with any semidefmite solver interfaced through 
YALMIP as follows: 

» solvesdp(F,h) ; 
» ysol = double (y); 
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5.12 Moment substitutions 



By performing explicit moment substitutions it is often possible to reduce significantly the 
number of variables and constraints in moment SDP problems. Version 2 of GloptiPoly 
implemented these substitutions for mixed-integer 0-1 problems only [3J. With version 3, 
these substitutions can be carried out in full generality 

GloptiPoly 3 carries out moment substitutions as soon as the left hand-side of a support 
or moment equality constraint consists of an isolated monic monomial. Otherwise, no 
substitution is achieved and the equality constraint is preserved. 

For example, consider the AW® Max-Cut problem studied in [3, §4.7], with variables 
Xi taking values —1 or +1 for i = 1, ... ,9. These integer constraints can be expressed 
algebraically as xj = 1. The following piece of code builds up the third relaxation of this 
problem: 

» W = diag(ones(8,l),l)+diag(ones(7,l),2)+diag([l 1] ,7)+diag(l ,8) ; 
» W = W+W; n = size(W,l); e = ones(l,n); Q = (diag(e*W)-W)/4; 
» mset clear 
» mpoK'x' , n) 

» P = msdp(max(x , *Q*x) , x.~2 == 1, 3) 
GloptiPoly 3.0 
Define moment SDP problem 
Valid objective function 

Number of support constraints = 9 including 9 substitutions 

Number of moment constraints = 
Measure #1 

Maximum degree = 2 

Number of variables = 9 

Number of moments = 5005 
Order of SDP relaxation = 3 
Mass of measure 1 set to one 
Total number of monomials = 5005 
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Perform moment substitutions 

Perform support substitutions 

Number of monomials after substitution = 465 

Generate moment and support constraints 

Generate moment SDP problem 



Moment SDP problem 

Measure label = 1 

Relaxation order = 3 

Decision variables = 465 



Semidefinite inequalities = 130x130 

We see that out of the 5005 moments (corresponding to all the monomials of 9 variables 
of degree up to 6), only 465 linearly independent moments appear in a reduced moment 
matrix of dimension 130. 

With the following syntax, moment substitutions are not carried out: 
» P = msdp(max(x'*Q*x) , x.~2-l ==0, 3) 

Mass of measure 1 set to one 

Total number of monomials = 5005 

Perform moment substitutions 

Number of monomials after substitution = 5004 

Generate moment and support constraints 

Generate moment SDP problem 



Moment SDP problem 

Measure label = 1 

Relaxation order = 3 

Decision variables = 5004 

Linear equalities = 6435 
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Semidef inite inequalities = 220x220 

Only the mass is substituted, and the remaining 5004 moments linked by 6435 linear 
equalities (many of which are redundant) now appear explicitly in a full-size moment 
matrix of dimension 220. 
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